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Two contemporary themes in the study of cold atomic gases are the creation of new 
exotic forms of quantum matter, and quantum simulations of systems already present in 
nature [USE]. By tuning the parameters for a collection of atoms and lasers one may 
address problems in quantum many-body systems or in high-energy physics |4j. In the 
first case degeneracy, quantum correlation, and entanglement are essential ingredients, 
whereas the latter case usually focuses on low-energy fluctuations of systems where a 
macroscopic fraction of particles reside in a single quantum state, often amenable to 
Landau descriptions. The versatility of Bose-Einstein condensates (BECs) allows the 
freedom to specify the geometry and topology of the order parameter to suit a particular 
purpose. For example, spinor BECs provide one way to realize order parameters with 
large symmetry groups, and hence exotic topologies 0 ei m ei ei dnum H 21 d3i nn tin] . 
It follows that the inclusion of spin-orbit coupling in such systems increases their 
complexity and introduces topological order H3 113 m ESI 1201 EH, a distinct 
classification for the order parameter. However, in order to access interesting physics 
and to simulate new regimes it may be necessary to extend beyond the usual notion of 
stability to metastable non-ground-state or non-equilibrium BECs. 

In this article, we develop some of the fundamentals underlying non-ground-state 
BECs in quasi-two-dimensional (quasi-2D) honeycomb lattices and the associated long- 
wavelength emergent theories [221 [23], EH ESJ 23 EH 29]. We focus in particular on 
superfluid fluctuations in the presence of Dirac points from a semiclassical perspective 
and by including lowest-order quantum effects. Quantum fluctuations are determined 
by solving the partial differential equations which describe dynamics of the low-energy 
modes for an arbitrary condensate profile. These equations are Lorentz invariant and 
comprise a relativistic generalization of the Bogoliubov-de Gennes equations (BdGE); 
thus we call them relativistic linear stability equations (RLSE). The RLSE provide a 
means of calculating vortex stabilities, yet their versatility extends beyond stability 
calculations to simulating a large number of established theories in addition to some 
exotic ones. This is because quasiparticles in BECs with inherent relativistic structure 
(e.g., linear dispersion, CVT invariance, multicomponent order parameter, etc.) can 
be tuned to have linear or quadratic dispersion with a zero or finite gap coupled to 
a condensate reservoir with a large number of possible internal symmetries.t In short, 
ease of construction and manipulation of BECs with quasi-relativistic dispersion and 
characteristically low sound speeds (on the order of centimeters per second) present 
ideal environments for simulating high energy phenomena. 

Moreover, the “no-node” theorem originally proposed by Feynman [36], which 
constrains conventional BECs, is circumvented for non-ground state (metastable) 
systems and in the case of spin-orbit coupling, as the order parameter in these systems is 
generally not positive-definite [3?, [38]. This property is a fundamental feature of quasi- 

f Relativistic effects occur in photonic systems as well. See for example Refs. [201 Ell S21E2]- Interesting 
dynamics in the merging of Dirac points for fermionic systems have also been investigated [341 135] 
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relativistic condensed matter systems. In particular, lifting the “no-node” theorem 
restriction leads to time-reversal symmetry breaking, which allows for exotic bosonic 
systems such as p-wave superfluids [39], chiral Bose liquids [40], complex unconventional 
BECs in high orbital bands of optical lattices @1] , and BECs with repulsive interactions 
that support bright solitons and vortices as well as skyrmions @211121E]. We point 
out that our system is identical to a quasi-2D BEC with spin-orbit coupling in either 
the long-wavelength limit or the strong tunable spin-orbit coupled limit, provided the 
interactions are also chosen to retain only the intra-component terms. To map to the 
strong spin-orbit coupled limit, however, the strength of the spin-orbit coupling term 
must be much larger than the quadratic term but still below the quantum critical point 
separating the spin-balanced and spin-polarized ground states m- 

Our results focus on three main topics: 1) physical parameters and constraints ; 2) 
linear stability of vortices] and 3) emergent theories. First, the physical parameters and 
necessary constraints to construct a non-ground-state condensate at Dirac points are 
explained in detail. The BEC is tightly confined in one direction and loosely confined 
in the other two directions. More precisely stated, magnetic trapping along the z- 
direction is such that excitations along this direction have much higher energy, by 
at least an order of magnitude, compared to the lowest excitations in the x and in¬ 
directions. Thus, an important step is to calculate the precise renormalization of all 
relevant physical parameters when transitioning from the standard 3D BEC to a quasi- 
2D system. In addition to this step we also account for renormalization due to the 
presence of the optical lattice potential which introduces an additional length scale 
from the lattice constant. We point out that microscopically the BEC obeys the three- 
dimensional nonlinear Schrodinger equation and we consider temperatures well below 
the BKT transition energy associated with two-dimensional systems. Nevertheless, 
throughout our work we often use “2D” for brevity, keeping in mind the quasi-2D picture. 
Condensation at Dirac points of the honeycomb lattice requires additional techniques 
beyond ordinary condensation, which we have detailed in our previous work |45j . In 
addition to the fields needed to construct the lattice one requires a resonant field which 
provides the time-dependent potential to “walk” atoms from the ground state (zero 
crystal momentum) to the Dirac point. The result is a transient configuration since a 
macroscopically occupied nonzero Bloch mode is not in thermodynamic equilibrium. 

Care must be taken when transferring atoms from the ground state to a Dirac point 
in order to minimize depletion out of the condensate. In general, one might expect some 
dissipation to occur due to secondary interactions within the condensate, and between 
condensed atoms and the lattice, quantum fluctuations, and thermal excitations, the 
latter two comprising the normal fluid. However, at the mean-field level repulsive 
atomic interactions within the condensate itself produce a single Hartree term which 
just shifts the total energy upward without causing additional depletion. Lattice effects 
are accounted for completely through the band dispersion, since we are not considering 
the presence of disorder or artificial impurities. Moreover, we consider only the zero- 
temperature case. There is certainly finite leakage into energetic modes lower as well as 
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higher than the condensate energy. However, such losses can only occur in the presence 
of higher-order dissipative terms in the Hamiltonian. In this article we restrict our 
analysis to the effects of first-order quantum corrections and apply our results to the 
special case of vortex background. 

The second major topic in this article addresses linear stability of vortices near 
a Dirac point. We first provide a detailed derivation of the RLSE then solve them 
for vortex solutions of the nonlinear Dirac equation |43j . The resulting eigenvalues 
determine the characteristic lifetimes of each vortex type. Solutions of the RLSE are 
inherently massless Dirac spinors with components that couple only through the Dirac 
kinetic terms. For a vortex background, RLSE solutions describe the quantum density 
and phase fluctuations near the vortex core. Physically, these are local undulations 
in the density profile, rigid translations of the vortex itself, and fluctuations in the 
speed of rotation. Although the latter is topologically protected, at the mean-field 
level quantum effects introduce small admixtures of different winding numbers into the 
vortex. These admixtures, which take the form of phase fluctuations, comprise the 
Nambu-Goldstone modes of the system. Near the vortex core they appear as bound 
states, the lowest of which are zero-energy modes (zero modes): static modes with zero 
energy associated with spatial translations of the center of the vortex. From a symmetry 
perspective, zero modes account for the fact that a vortex breaks the translational and 
rotational symmetry of an otherwise uniform system. We will address the various modes 
in generality when we discuss the associated reductions of the RLSE to other well known 
equations. 

Our work culminates with the connection to several other important areas of 
physics including relativistic Bardeen-Cooper-Schrieffer (BCS) theory. In addition to 
continuous space-time dependence, quasiparticle solutions of the RLSE are labeled by 
two indices associated with the lattice pseudospin valley and the particle-hole structure 
analogous to Nambu space from BCS theory. In order to avoid confusion we will refer 
to these as valley and Nambu indices, and reserve the particle-hole terminology to 
distinguish between the two states in either the valley or Nambu space. The RLSE are 
formulated to describe excitations in a repulsive Bose gas but can be reinterpreted 
as excitations in a theory comprised of attractive particles upon pseudospin valley 
particle-hole exchange. This symmetry is a consequence of the combined symmetry 
of charge conjugation (C), parity transformation (V), and time reversal (T), which is 
fundamentally related to the structure of the Dirac operator. Retaining a mass gap, 
an intermediate step in Dirac-point condensation 051. and adding atomic interactions 
between nearest-neighbor lattice sites |46j extends the RLSE to the Dirac-Bogoliubov- 
de Gennes equations (DBdGE) [4?; 08, 09j [50 ]. provided valley particles and holes 
are interchanged. This connection is significant as DBdGE are required for a broader 
description of superconductivity beyond the standard BCS formalism, particularly for 
superconductors with a high Fermi velocity. Indeed, a relativistic formulation of BCS 
becomes important for elements with large atomic number [Z, > 40), in neutron stars 
where superfluidity is expected to play a major role in “glitches” mm [53J . and in 
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Figure 1 . (color online) Schematic overview for applications of the RLSE. The RLSE 
are derived for either a metastable BEC at Dirac points of a honeycomb optical or for 
a strong spin-orbit coupled BEC. Various reductions of the RLSE or augmentation by 
the addition of a mass gap and nearest-neighbor interactions (lattice case), or general 
spinor couplings (spin-orbit coupled case), leads to five categories of equations from 
different areas of physics. 


color superconductivity where the strong nuclear force provides the attraction between 
fermions In the nonrelativistic and semiclassical limits the RLSE reduce to the 
Andreev equations. These equations were originally formulated to address physics 
of nonuniform superconductors, for instance a type-1 superconductor near a normal- 
superfluid interface or a vortex in a type-II superconductor [551156] . Interestingly, we find 
that the RLSE reduce to the Majorana equation inside the core of NLDE vortices. From 
a fundamental standpoint the Majorana equation describes relativistic fermions that are 
their own anti-particle ra- In condensed matter systems, and in particular our problem, 
finite-energy phase fluctuations inside the core of an NLDE vortex connect smoothly 
to Majorana zero modes. This is significant as Majorana zero modes are presently 
of great interest in such fields as topological quantum computation [58], topological 
insulators [59], and more generally in the study of non-Abelian anyons and fractional 
statistics mmmm- Figure [I] provides a schematic overview of some of the theories 
and physical regimes encapsulated in the RLSE. 

This article is organized as follows. In Sec. [2] we discuss physical parameters, 
constraints, and regimes. In Sec. [3] we analyze superfluid excitations for a Bose gas in the 
honeycomb lattice from a semiclassical perspective. Section [4] contains two derivations of 
the RLSE according to paths dictated by two possible orderings of the tight-binding and 
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continuum limits of lattice Bloch functions. In Sec. [5J stability analysis is performed for 
vortex solutions of the nonlinear Dirac equation by solving the RLSE for quasiparticle 
functions and eigenvalues. In Sec. [6j we examine several reductions of the RLSE to 
other well known equations. We map the RLSE to the equations for relativistic BCS 
theory and demonstrate the non-relativistic limit to standard BCS theory. In Sec. [7J we 
conclude. 

2. Renormalized Parameters and Physical Constraints 

To obtain the correct renormalized parameters for the NLDE we proceed by two steps. 
First, we follow the transformation of the 3D NLSE parameters as we reduce to the 2D 
NLSE. Second, we take the long-wavelength limit of the 2D theory at the Dirac point 
to get the NLDE, which induces a second renormalization of the parameters. 

2.1. Transition from 3D to 2D nonlinear Schrodinger equation 

A BEC comprised of N atoms of mass M is described by a wavefunction ^(r,t) 
which solves the time-dependent nonlinear Schrodinger equation. The single-particle 
density is defined as |-0(r,t)| 2 , the BEC density p(r,t) 2 = A/j^(r,t)| 2 , and the phase is 
(f) = arg[-0(r, t)], with the superfluid velocity given by v s = K\7f>/M. The two-particle 
interaction strength is g = Anh 2 a s /M and the healing length is £ = 1/\/8nha s , where a s 
is the s-wave scattering length for binary collisions between atoms. We take a s > 0 so 
that g > 0, i.e., we consider only repulsive interactions, leaving attractive interactions 
for future studies. Throughout our work, we treat the case of an axisymmetric system 
associated with a harmonic trapping potential with two large dimensions described by 
a radius R = \J x 2 + y 2 , and a small dimension transverse to the plane described by 
the length L z . The average density which appears in £ is then n = N/(ttR 2 L z ). Note 
that ip(r,t) has dimensions of length -3 / 2 so that g has dimensions of energyx length 3 . 
Another important quantity is the speed of sound in the condensate, which is defined 
as c s = y/gn/M. 

Transforming to the 2D regime requires that a, « < ( |3, 65], which ensures 

that the condensate remains in the ground state in the transverse direction, and L z <C R , 
which ensures that excitations along the plane have much lower energy than those in 
the transverse direction. The wavefunction can then be separated into longitudinal and 
transverse modes, following similar arguments as in Ref. (2a 

V>(r, t) = (AL,y 1/2 f(x, y)h(z)e ~ v ‘ ,h , (1) 

where f(x,y) and h(z) are the dimensionless spatial functions that describe the 
longitudinal and transverse normal modes, respectively, and p is the chemical potential. 
Projecting onto the ground state of the transverse dimension h gs (z), gives us an 
effectively 2D wave equation. In the case where L z ~ £, h gs (z ) is just the 
ground state of the one-dimensional particle-in-a-box solution [22], we then have 
h gs (z) = sin(7T z / L z ). It may be convenient to express L z and R in terms of 
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the trap frequencies Ld X: Lu y , and u z , in which case we may write L z = (h/MuJz ) 1 ^ 2 , 
R = -^/hM~ l {l/u x + 1/c Oy). The transformation is then completed by defining the 
renormalized 2D chemical potential and interaction as 


/Dd = 9 + 


h 2 7T 2 

2 ML 2 Z ’ 


92 D 


3^ 

2 L z ‘ 


( 2 ) 


Through this process, we obtain a reduction of the 3D nonlinear Schrodinger equation 
to a 2D form with renormalized chemical potential and interaction given by Eq. (J2]). 
The 2D renormalized average density can be related to the 3D average density using 
the transverse oscillator length or frequency 

1/2 

n. (3) 


N 


n 2D = ~ Lz n ~ 


h 


Mix , 


Using this definition and the 2D single-particle wavefunction, if(x, y ) = A~ 1 ^ 2 f(x, y ), we 
can write the 2D condensate density as p 2 n(x,y) = N |'0(:r,y)| 2 . The 2D renormalized 
healing length can also be constructed which we find acquires only an extra numerical 
factor 


2\ 1/2 1 


^ 2D 3 J y/8irna s 



(4) 


Similarly, we find the 2D speed of sound to be c S 2 D = \JgxD^n/M = (3/2)^ 2 c s . 
It is important to keep track of the effect of the reduced dimensionality on the 
dimensions of the constants: //>( x , y) now has dimensions of length -1 , (/ 2 D has dimensions 
energyxlength 2 , and n 2 B has dimensions length -2 . 


2.2. Derivation of nonlinear Dirac equation from 2D nonlinear Schrodinger equation. 


The derivation of the nonlinear Dirac equation begins with the second quantized 
Hamiltonian for a 2D system with the bosonic field operators if = = ^{x^yfl) 

obeying bosonic commutation relations in the Heisenberg picture. We then expand in 
terms of Bloch states belonging to A or B sites of the honeycomb lattice which breaks up 
the bosonic field operator into a sum over the two sublattices. The spatial dependence in 
this expansion is encapsulated in the exponential Bloch wave and the Wannier functions 
w(x,y) which are then integrated out leaving only number-operator terms in the form 
of a Dirac-Hubbard Hamiltonian (a detailed derivation is presented in Ref. 122) 


h = -t h £ 

<A,B> 


a r be ik<TA - TB) + tfae- ik{rA ~ rB) 


+ y ^ a^d^aa + y ^ tf&bb. (5) 

A B 


The bracketed A and B summation index signifies a sum over nearest-neighbor A and B 
sites, accounting for inter-sublattice hopping through the individual sublattice creation 
and destruction operators af and a,b, respectively, with an accumulated phase which 
depends on the separation vectors r ab = — r^. The interaction terms appear quartic 

in operators of the same site, as contact interactions are local and do not couple the 
A and B sublattice. Note that th and U are the hopping energy and 2D renormalized 
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interaction, respectively. Terms proportional to ad a and Db just count the total number 
of atoms in the system, and have been neglected in Eq. (J5| as an overall constant. 
Equation ([5]) is the Hubbard Hamiltonian divided into two degenerate sublattices A 
and B, appropriate to the honeycomb optical lattice. 

The time evolution of the bosonic operators a and b is computed according to the 
standard Heisenberg prescription: ihd t a = [a, H], Although we suppress the lattice 
site indices, implicitly the operator d hJ which destroys a boson at site (i, k ) satisfies the 
bosonic commutation relation [d;j, a\, ■,] = 5a>8jji. Using this fact and a similar relation 
for b iJ} we compute the commutator of a with the Hamiltonian Eq. ([H]) to obtain 


ik(r Ai .-r Bi .) ^ l r,. r r fl( . 3) _ ni ) £ ^ p ^(r A . . -r B( . ) 


ihdthij = - t h bije + b ^ j) _ ni < 

T Ud t j 


n 2 


( 6 ) 


where the first three terms on the right hand side represent transitions from the three 
B-sites nearest the ( i,j) th site of the A sublattice and ni and ri 2 are primitive cell 
translation vectors for the reciprocal lattice, as shown in Fig. [7]^a) . The same method 
yields a similar equation for the B sublattice. We next calculate the time rate of change 
of the expectation value of Eq. ([TTJ) with respect to on-site coherent states. For the 
(i,j) th site this is defined as 



n 


which is an eigenstate of the destruction operator d lJ with eigenvalue a lyJ e C, as can 
readily be verified. Taking the full wavefunction to be the direct product over all MxM 
lattice sites (including the B sublattice), we have 


MxM 


I a, b ; M x M) = (^) (ja^) (g) . 


( 8 ) 




From the definition of the single-site coherent state Eq. (J7|), one may then verify that 


djj|a, 6; M x M) = a it j\a, b\ M x M ), 


(9) 


and 


a, b ; M x M\a\ • = (a, 6 ; M x M\a* ■ 




( 10 ) 


Inserting the nearest-neighbor vectors 5i, S 2 , and ^3 into the exponentials in Eq. (11) 
Fig. 0 (b)), and doing the same in the equation for b, then taking the expectation 


see 


with respect to the state Eq. (| 8 j) we obtain coupled equations of motion for discrete, 
on-site, complex-valued amplitudes 


iha id = -t h [bi,je lk ' S3 + %j)_ ni e * k ' <51 + %j)-n 2 e* k ' <52 ] + Ua* tj a itj a id , 
ihb id = —t h \a itj e~ l ^ h + a^ j)+ni e~ lk ' Sl + a^ j)+n2 e _lk ' <52 ] + Ub^biJ)^ 


( 11 ) 

( 12 ) 


The NLDE is derived around the linear band crossings between the A and B sublattices 
at the Brillouin zone corners [ 66 ], called the Dirac cones in the graphene literature ra. 
To this end, we insert appropriate values for the the nearest-neighbor displacement 







Superfluid fluctuations and emergent theories 


9 


vectors 5 and evaluate the wave vector k at the Brillouin zone corner, defined by 
k = K = (0,4tt/ 3) , ^ = (l/2x/3, — 1 / 2 ), <J 2 = (1/2^, 1/2), S 3 = (-1/^3, 0). Finally, 
taking the long-wavelength limit of the resulting equations recovers a continuum theory 
but with a Weyl spinor wavefunction 4/ = {ip a-, J’b)- 

The key point in discerning the correct normalization (and thus other related 
quantities) is the contraction of the many-body bosonic operators between localized 
coherent states. The parameters \a u j ( 2 and | 6 j . ; | 2 which label the coherent states for 
the A and B sublattices at site {i,j), emerge as the number of atoms at each site, 
so that a i} j and b t ,j become continuous amplitudes ip A { r ,t) and ipB{ r ,t ) in the long- 
wavelength limit, as stated above. Note that the complex moduli of these amplitudes 
are pure dimensionless particle numbers, not densities, since they result from taking the 
spatial integral over the lattice. With the area per lattice site given by Ai = \/3a 2 /4, 
the local time-dependent sublattice densities can be reconstructed as: p A (B){*,t) = 
\ipA(B){ r fl)\ 2 /Ai. Then, the dimensionally correct sublattice mean-held wavefunctions 
must be given by ip A (B){^, t)/y/Ai = (16/3a 4 ) 4 / 4 ip A ( B )(r, t), where a is the usual lattice 
spacing. The correct normalization procedure can now be deduced by writing down the 
total number of particles in the system 

/*27T 

N = (16/3C1 4 ) 1 / 2 / dtp 
Jo 

where the upper limit of the radial integral is taken large enough so that the integrand 
is negligible. The total number of atoms of the system, N, appears on the left-hand 
side. 


Jrdr{\ip A {r,<j)-,t)\ 2 + \ip B {r,<f>;t )\ 2 ), (13) 


The 3D to quasi-2D reduction and continuum regime result in an effective atomic 
interaction U, a renormalized version of the usual interaction g. We arrive at the 
explicit form for U by first approximating the lowest band on-site Wannier functions by 
the ground state of the harmonic oscillator potential. Integrating over the area of one 
site, we obtain a new local interaction strength 

2 


U = #2D 


\/3a 2 


n 2Y) dxdy\wi{x,y)\ 


— 92 D 


V?>a 2 


— 2 


n; 


2D 


27 t£ 2 J 


(14) 


where £ is the oscillator length of a lattice potential well. It is often more practical 
to express the area of one site in terms of the lattice constant ni 2 = v / 3o 2 /4, and all 
other parameters in terms of the corresponding 3D parameters. Using Eqs. (§-<§ , the 
interaction takes the form 

U = L z g n 2 a • (15) 

Note that U has dimensions of energy. 

We can now identify the main parameters which appear in the NLDE. The 
dimensionful coefficient which multiplies the Dirac kinetic term is the effective speed of 
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Parameter 

Symbol/Definition 

Value 

Range 

Plank’s constant 

h 

1.06 x 10 34 j ■ s 


Boltzman’s constant 

ks 

1.38 x 10 23 j ■ KX 1 


Mass of S7 Rb 

M 

1.44 x 10 ~ 25 kg 


Number of atoms 

N 

3.00 x 10 4 

o 

1 

o 

o 

Wave number of laser light 

k L 

7.57 x lO 6 !!^ 1 

4.19 x 10 6 -4.19 x lO’m - 1 

Lattice constant 

a = 47 r/ 3 fci 

0.55 gm 

0.30 — 0.70 /im 

Recoil energy 

E r = h 2 k 2 L /2M 

0.16 gK 

0.049 - 4.90 

Lattice potential 

Vo — 16-Ei? 

10.1 gK 

0.79- 10.1/iK 

Hopping energy 

4 = 1-861 (V 0 /E R ) 3/i E R e- 1 - 5S2 v' Vo/El ‘ 

16.8 nK 

3.49 nK - 1.90 ^iK 

Scattering length 

CL S 

5.77 nm 

5.00- 10.0 nm 

Average particle density 

n 

5.86 x 10 18 m - 3 

10 15 - 10 21 m - 3 

Two-body interaction 

g — Airh 2 a s /M 

41.0 K ■ nm 3 

22.36- 52.18 K ■ nm 3 

Healing length 

£ = l/y/8ir na s 

1.10 /zm 

< 1.50 ^m 

Sound speed 

c» = pgfi/M 

4.82 x lO -2 cm/s 

5.83 x HT 3 - 0.825 cm/s 

Sound speed (2D) 

c« 2 D = (3/2) 1/2 c s 

5.90 x lO -2 cm/s 

7.14 x HT 3 - 1.01 cm/s 

Healing length (2D) 

6 d = (2/.3) 1/2 £ 

1.75 /im 

< 2.45 (im 

Transverse trap energy 


22.17 nK 

0.21 - 56.5 nK 

Transverse oscillator length 

L z = (h/Mui z y/ 2 

1.50 fin l 

< 3.0 41 m 

Average particle density (2D) 

n 2 D = L z n 

4.50 x 10 12 m - 2 

10 9 - 5.00 x 10 15 m -2 

Effective speed of light 

ci = thCtVs/zh 

5.31 x 10 ~ 2 cm/s 

< 5.40 x lO -2 cm/s 

Dirac kinetic coefficient 

ci = hci 

2.07 nK • /im 

< 5.72 nK ■ 41 m 

Dirac nonlinearity 

U — L z gn 2 3\/3 a 2 / 8 

1.07 nK 

< 2.36 nK 

Dirac healing length 

^Dirac ~ t/ l fl\/3/2t/ 

3.80 fim 

0.50 - 50 . 04 im 


Table 1 . Physical Parameters for the NLDE typical for a BEC of 87 Rb atoms. The 
renormalized parameters are expressed in terms of fundamental quantities. The range 
of possible values account for the physical constraints discussed in the main text. 


light Ci ~ 5.31 x 10 -2 cm/s (compare to the analogous coefficient for relativistic electrons 
c ~ 3.00 x 10 8 m/s). In terms of fundamental constants we find q = thaV3/2h, where 
a is the lattice constant and th is the hopping energy. The natural length scale of the 
NLDE is the Dirac healing length £oirac = hci/U = tha\/3/2U , which characterizes the 
distance over which a disturbance of the condensate will return to its uniform value. 
We see that £oirac has the correct dimension of length. To simplify the notation, for the 
remainder of our paper we will omit the 2D subscript on all parameters with typical 
values as can be achieved in present experiments |T5j- Finally, the quantity U which 
appears in the NLDE determines the strength of the nonlinearity. We have provided a 
full list of relevant parameters associated with the NLDE in Table [l] 

2.3. Physical constraints 

The realization of the NLDE in a condensate of 8 'Rb atoms requires that several 
constraints are satisfied which we now list and discuss: 

(i) Landau Criterion. In order to avoid the instabilities associated with propagation 
faster than the sound speed in the condensate, we require that the effective speed 
of light is less than the 2D renormalized speed of sound. 

(ii) Long-wavelength Limit. The NLDE describes propagation of the long-wavelength 
Bloch envelope of a BEC near the Dirac point. Thus, a necessary condition for 
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realizing the NLDE in the laboratory is that the Dirac healing length must be 
much larger than the lattice constant. 

(iii) Relative Lengths for 2D Theory. In order to obtain an effectively 2D system, the 
vertical oscillator length must be much smaller than the trap size along the direction 
of the plane of the condensate. 

(iv) Relative Energies for 2D Theory. Analogous to the previous restriction, this 
condition relates to the 2D structure but pertains to the energies of the system. The 
key point is that we must avoid excitations vertical to the plane of the condensate 
while enabling them along the plane: the chemical potential and temperature must 
be less than the lowest transverse excitation energy. 

(v) Weakly Interacting Regime. The NLDE and RLSE are derived for a weakly 
interacting Bose gas. This ensures both the stability of the condensate as well as 
the effective nonlinear Dirac mean-field description. We then require the interaction 
energy to be significantly less than the total energy of the system. 

(vi) Dirac Cone Approximation. For a condensate in the regime where the NLDE 
description is valid, we require that the linear approximation to the exact dispersion 
remain valid. As in the case of graphene, large deviations from the Dirac point 
induce second order curvature corrections to the dispersion. Thus, we must quantify 
the parameter restrictions which allow for a quasi-relativistic interpretation.t 

(vii) Lowest Band Approximation. We derive the NLDE and RLSE assuming that the 
lowest band is the main contribution to the dispersion. 

Having stated each constraint, we can now address each one in detail and explore the 
conditions under which each is satisfied. In the following, we consider a BEC comprised 
of 8 'Rb atoms where all numbers used are listed in Table [l] and are experimentally 
realistic ra- First, the Landau criterion pertains to the effective velocities in the 
BEC. Stated mathematically, the Landau criterion requires that q/c s 2 d < 1- Using 
the definitions for the effective speed of light and the sound speed found in the first part 
of this section, we compute q/c s 2 d = 0.90, which satisfies the inequality. 

The length constraints are as follows. The long-wavelength limit is defined by 
^Dirac/o 1, for which we find that £oirac/a = 6.91. For an effectively 2D system, the 
required length constraint implies the condition L z <C R. Taking R « 100 a (a typical 
condensate size), and using a realistic value for the vertical oscillator length (Table [I]), 
we obtain L z = 2.73 a, which satisfies the constraint. Moreover, we require a healing 
length close to or less than the transverse oscillator length. With £ = 1.10 pm and 
L z = 1.50 pm, we find that this condition holds. 

The energy constraints may be stated as p, ksT <C hui z . We can solve the 
NLDE for the lowest excitation to obtain an expression for the chemical potential 
p = hcik + U |T| 2 [TlJ. Next, we evaluate this expression using the lowest excitation 

f We note that the Dirac cone approximation is not necessarily adhered to in analogous honeycomb 
photonic lattice systems. See for example Refs. [U51 ITjU] . 
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in a planar condensate of radius R ~ 100a, which has wavenumber k ~ ir/2R = 
2.86 x 10 4 m _1 . The interaction U is computed using Eq. (15) for the binary interaction 
g and mass M pertaining to a condensate of s 'Rb atoms. Finally, for a uniform 
condensate we take |T | 2 « 4/\/3 (Eq. (13)) and the constraint on the chemical potential 
becomes p = 2.59 nK < 22.17nK, which is satisfied. For the temperature, we require 
T <C hjj z /kB- Using the data in Table [l] for the vertical oscillator frequency, we obtain 
the upper bound for the temperature T <C 22.17 nK. This is a reasonable requirement 
given that BEC occurs for T in tens or hundreds of nanoKelvins or as low as picoKelvins. 

Next, we examine constraints on the particle interaction. To check that we are in 
the weakly interacting regime, i.e., that U / p -C 1 , we use the value for the chemical 
potential p which we have just computed and compare this to the interaction energy 
U, whereby we find that U / p = 0.41. An essential feature of the NLDE is that 
characteristic fluctuations are close enough to the Dirac point so that the linear Dirac 
cone approximation remains valid. Expanding the exact dispersion near the Dirac point, 
we obtain p{k) = U ± th {aV?>k/2 + a 2 k 2 /8 — a 3 -\/3fc 3 /48 + ...), where k is a small 
deviation away from the Dirac point. The first term gives the linear Dirac dispersion. 
Higher order corrections describe curvature of the band structure away from the Dirac 
point. From the second order term we see that the NLDE description is valid for 
ak/\J8 -C 1. This determines a lower bound on the wavelength for fluctuations of the 
condensate: A m i n 3> (27r/v / 8)a. Linear dispersion places an additional constraint on 
the chemical potential: \p\ <C U + 6 1^ ~ 101.9nK. From the value of the chemical 
potential already obtained, we find p = 2.59 nK <C 101.9 nK. Finally, weak short range 
interactions at very low temperatures justifies a lowest-band approximation to describe 
the physics of the NLDE. 


3. Superfluid excitations near a Dirac point 

The mean-field physics of single-particle states for a collection of fermions with Fermi 
energy near a Dirac point of a honeycomb lattice has been studied exhaustively 
and is discussed in various comprehensive articles UM EH Eg. For systems of 
bosons, however, one must carefully consider the meaning of condensation in the 
presence of Dirac points. To discuss BECs and Dirac points together one must address 
the compatibility of single-valuedness for phase functions required for stable vortex 
formation in a proper superfluid description, with the half-angle phase winding when 
circumnavigating a single Dirac point, i.e., the geometric or Berry phase m 

3.1. Geometric and dynamical phase structure 

To address these issues, we first review some relevant information treated in most review 
articles on graphene, as this information is true for cold bosonic atoms as well [27]. The 
single-particle spectrum of the honeycomb lattice exhibits zero-points, or Dirac points, 
in the reciprocal lattice associated with crystal momentum K = (0, ±47r/3a) rotated 
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Figure 2. (color online) Bragg wave structure at the Dirac point of the honeycomb 
lattice. Constructive interference of Bragg reflected waves for the Dirac point 
wavevector K (see Fig.[7](c)) produces wavefront density peaks at A (red) and B (blue) 
sublattices. 


by 0,27r/3, 47 t/ 3 , where a is the lattice constant shown in Fig. [2j Dirac points occur 
when the crystal momentum is tuned to the natural periodicity of the lattice with 
standing waves established due to Bragg scattering of the wave function. Reflection at 
the Brillouin zone edge is shown in Fig.[2j where one adds up projections of the vectors rix 
and n 3 along the direction of the crystal momentum vector K connecting points on the 
A sublattice of equal phase, Ai and A 2 , to a third point A 3 . In particular, at the Dirac 
point this sum results in a net 2tt accumulated phase angle at A 3 . I 11 Fig. [2j the A and B 
sublattice wavefront density peaks are shown as red and blue dashed lines, respectively. 
In the tight-binding limit, the full lattice Hamiltonian reduces to two operators which 
couple the degenerate triangular A and B sublattices. The single-particle dispersion is 
computed by solving the 2 x 2 eigenvalue problem in momentum space determined by 
the Hamiltonian 


H(k) 


0 E (k)* \ 
E( k) 0 ) 


(16) 


where the matrix elements come from computing the sublattice hopping energies 

E( k) = -t h (1 + e ini k + e in2 k ) = \E{k)\e~ iW . (17) 


Specifically, one finds 


^1 + 4 cos 




+ 4 cos 2 



|£(k)| = t h 


(18) 
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The eigenfunctions of H in Eq. (|16|) are 

*±(k) = 


e i£0(k) 
e i(-€+i)0( k ) 


(19) 


with eigenvalues ±|£'(k)|. Physically, the parameter i e [R comes from an extra 
U(l) phase degeneracy and reflects the gapless symmetry of the system under spatial 


translations of the atomic density at the Dirac point. The matrix in Eq. (16) describes 


the amplitude and phase associated with real-particle tunneling between neighboring 
lattice sites. In particular, the phase of the wavefunction gets multiplied alternately by 
factors of e ±l ^\ so that no net phase is accrued when circumnavigating a closed path in 
the lattice. In contrast, long wavelength modes propagating in the lattice are described 
by linearizing the phase angle 0(k) so that the local lattice scale variations in the phase 
structure are neglected, in which case one should expect a net phase accumulation. 

We are particularly interested in this net geometric or Berry phase since we must 
factor it into the phase winding for vortex solutions of the NLDE. Although most 
treatments of the subject use a momentum space argument, here we use instead a more 
direct analysis in real space. We expand the Hamiltonian and eigenstates near the Dirac 
point by taking k = K + hk, with K = (47r/3a)y and <5k the small expansion parameter, 
i.e., we consider small deviations from the sublattice Brillouin zone corner. In real space, 
this amounts to a derivative expansion of Eq. 0 in terms of the directional derivatives 
ri! • V and n 2 • V. The first-order term gives the massless Dirac Hamiltonian and Dirac 
equation while higher products of derivatives provide corrections that probe the 
finer details of Bragg scattering around the Dirac point. To isolate the geometric phase, 
we consider adiabatic transport around a closed loop. Adiabaticity ensures that we 
do not accumulate a dynamical contribution to the phase and restricts the path to 
energy eigenstates states nearest the Dirac point. A direct way to accomplish this is by 


linearizing Eq. (16) in real space, solving for the eigenstates in plane-polar coordinates 


r and 9, and restricting to paths with large radii R = r a. The large radius limit 
allows us to access only the longest wavelength modes that vary mainly tangentially 


with minimal radial contribution. Equation (|16|) then reduces to 

—e~ l6 d e 

0 


H(6) 


ihcj 


0 


R 


f e do 


from which we find the eigenstates 


*±{ 0 ) = 


e-ie/2 

± e i6 ' 2 


( 20 ) 


( 21 ) 


and energies ±ftw/2, with u = cjR and q the effective speed of light. Note that in 


this limit the degeneracy in Eq. (19) is lifted and eigenstates are forced into the form 


Eq. (21), which acquires a net phase of tt under a full 27 t rotation. Thus, linearization 


of Eq. (16) leads to a double wrapping of the phase angle 0 around the polar angle 9. 


In the case of a vortex, we can compensate for the Berry phase by requiring 
half-winding in the overall dynamical phase multiplying the spinor order parameter: 
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exp [i(£ + 1/2)0], As a result, the geometric phase becomes identified with the relative 
phase between the two sublattices. Hence, stable vortices are required to have half¬ 
integer internal geometric winding plus an overall half-integer dynamical winding such 
that the superfluid velocity is the sum of the gradients of both phases, effectively splicing 
together the internal and external phase. 


3.2. Superfluid regime and dimensional analysis 

To address superfluidity in the honeycomb lattice using a semiclassical approach, 
consider a thermal excitation with crystal momentum p (measured from the Dirac 
point) interacting with an atomic gas at the Dirac point (poirac = 0) producing an 
excitation in the gas with momentum p'. It follows from energy conservation, A E = 0, 
that 


A E = ± 


IP ~ P 

2m* 


Pi 


=FQ|P-P'I +-E(p') =F ; ± Ci|p| = 0. 


( 22 ) 


In Eq. (22), for generality we assume that the incoming thermal mode may be in a Bloch 


state far enough removed from the Dirac point such that second-order corrections are 
important. Thus, m* is the effective mass related to the dispersion curvature, E( p') 
is the energy of the quasiparticle excitation in the gas, and the upper and lower signs 
refer to negative and positive dispersion branches, respectively. We first examine the 
linear regime for which p zs p' « cim *, in which case we can neglect quadratic terms 
in Eq. (22). Keeping only linear terms and using | p — p' | = \J p 1 — 2pp'cos9 + p' 2 , with 


p — |p | and 9 the angle between p and p', Eq. (22) forces the constraint 

E{ p') 


COS0 = ± 


Cip' 


(23) 


and four conditions determined by the different sign combinations for incoming and 
scattered modes. When the signs of the incoming thermal and scattered condensate 
modes are the same we find that 0 = 0. On the other hand, if the energies of incoming 
and scattered modes have opposite sign, we obtain 0 = ir, thus scattering in the reverse 
direction occurs between Dirac particles and anti-particles as one should expect. Notice 
that in this linear regime conservation of energy places no additional constraints on 
p and p', so that in our mean-held analysis an equilibrium between condensate and 
non condensate atoms is maintained: the incoming mode transfers all of its energy and 
momentum to an excitation of the condensate leaving a single outgoing excitation with 
the same energy and momentum. 

A second regime of Eq. (22) corresponds to the condition p' << p < cpm*, which 
yields the constraint 


p cos 0 = ± m 


M pO 

pi 


(24) 


leading to the same conditions as in the previous case for the scattering angle 0, but now 
with an additional constraint for an upper bound critical momentum p c below which no 
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excitation can be created in the condensate 

W)l 


p c = m Him p / - 


p' 


(25) 


Equation (25) recovers Landau’s criterion for superfluidity but here in terms of the 
absolute value of the quasiparticle energy to account for scattering into negative energy 
states. We point out that the absolute value in Eq. (25) is a strictly a consequence 
of energy conservation and the presence of quadratic terms in Eq. (22); the various 
sign combinations in Eq. (24) are taken into account through the scattering angle 6. 
With E( p') = ±cjp, the upper critical bound is just p c = m*ci. Below this value (and 
for p » p') thermal modes cannot interact with the condensate, thus superfluidity is 
preserved. For p > p', however, as expected we see a breakdown of superfluity. In our 
analysis we have nowhere included details of the interaction; only a knowledge of states 
near the Dirac point was needed. Once we consider quantum effects and details of the 
interaction our results will change significantly, as we will see in Sec. [5j 

Consider again a non-condensate excitation with initial momentum p interacting 
with the condensate by transferring all of its momentum and energy to the condensate. 
A secondary excitation is then emitted with exactly the same momentum and energy. 
Since the initial and final excitations are indistinguishable, we can view this process as 
a single excitation interacting weakly with the condensate and continuing on its way 
with only an average self-energy correction. At the quantum level and to first order in 
the interaction U, a single interaction point must be averaged over the volume (area in 
quasi-2D) of the condensate. Since we are dealing with very long wavelengths the result 
is a nonlocal collective excitation formed as a composite of the initial incoming mode 
dressed by the condensate background. At long wavelengths linear perturbation couples 
particles and holes, which amounts to reducing the power in the dispersion relation E(k) 
by a factor of 1/2. In contrast, at shorter wavelengths (higher energies) the incoming 
excitation couples with a quasiparticle locally, so that the available states for thermal 
and condensate modes remain distinct. Elaborating further, the number of accessible 
states less than k for the undressed excitation plus condensate is f l(k) = a r k r , and the 
dispersion is E = ±cihk s . Here for our argument we leave the constants r, s, and a r 
general. Thus, Q(E) = 2a r E r I s /(hci) r ! s where the extra factor of 2 accounts for both 
positive and negative eigenvalues. This yields the density of states 

d r /AD 5 ) - 1 

< 26 > 

In order to maintain D(E ) constant when transitioning from short to long wavelengths, 
Cip > U —y Cip < U, imposing the 1/2-power reduction in the dispersion, i.e., 
fl(fc) = a r k r —> a r k r / 2 , requires that we also take s —> s/2. The renormalized energy is 
then E oc ifc 1 / 2 (for s = 1). The proportionality constant must involve the quasi-2D 
interaction U which we determine through dimensional analysis to be 

E(p) = ±^/U^p. (27) 

Note that Eq. (27) leaves out the possible form E(p) = ±y/—Ucip, which displays 
a low-momentum dynamical instability. However, this is regularized by accounting for 
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a finite system size which imposes a lower momentum cutoff |p| m m = 27 tH/R, where R 
is the radial size of the condensate. For the usual harmonic trap with frequency hi we 
have R = (h/ . By dimensional analysis, in terms of the quasi-2D renormalized 
average particle density n and interaction U, we obtain the stability requirement for 
the oscillator length R < y/ci‘nh/{nU). From a practical standpoint the lower bound 
|p| min removes the longest wavelength modes which opens an insulating buffer between 
the positive and negative parts of the spectrum in addition to regulating the dynamical 
instability. 

4. Relativistic linear stability equations 

Bogoliubov’s method was originally introduced in his 1947 paper ca (see also CHI EH] for 
thorough contemporary treatments), and the concept later generalized by Fetter [80] to 
accommodate nonuniform condensate profiles. The latter formulation gives a convenient 
method for computing quasiparticle states and the associated eigenvalues by substituting 
the spatial functions for a particular background condensate into a pair of coupled 
differential equations, and then solving the resulting eigenvalue problem. Fetter’s 
extended method was designed with a vortex profile in mind, and has proven successful 
for computing stability of vortices in trapped condensates, but also for gaining a deeper 
understanding of general vortex dynamics [EH E2J E3]. The set of equations that 
we derive in this section form the counterpart to Fetter’s equations, but for trapped 
condensates that exhibit Dirac points in their dispersion uu. We call them relativistic 
linear stability equations because of the quasi-relativistic context here and the similarity 
to equations that appear in relativistic fluid dynamics. It is noteworthy that our result 
is not limited to the honeycomb optical lattice but applies generic-ally to any system 
where the linear dispersion and bipartite structure are present, and where the contact 
interaction between constituent bosons is weak. 

Our derivation of the RLSE relies fundamentally on Bogoliubov’s method as 
the underlying principle, and refers to Fetter’s work [811] for technical considerations 
regarding nonuniform condensates. First, we recall the second-quantized many-body 
Hamiltonian for weakly interacting bosons 

H = j d 2 r $ H 0 ip + | J d 2 r $ $ , (28) 

(29) 

where 

H 0 = -^-V 2 + V( r). (30) 

2m 

Here, V (r) is the lattice potential and g is the strength of the contact interaction. 
The first step is to decompose the wavefunction as the sum -0(r) = C( r )®o + </>(r), 
where we have split the wavefunction into a part that describes the condensate (first 
term) and satisfies the bosonic commutation relation [a 0 , <2 q] = 1, and a second 
part that describes small quasiparticle fluctuations. The operator in the first term 
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destroys a particle in the mean-field (, which, by itself, is a good approximation to 
ip. The second term destroys a particle in a number of single particle basis states 
of the noninteracting system, and describes the part of ip that deviates from the 
mean field. Taking the Bogoliubov limit requires do —> N 0 , where No is the total 
number of condensed atoms, but we choose to compute the commutator before taking 
this limit in order to retain the effect of the presence of a macroscopic condensate 
field. We can obtain the commutation relations for <p and by knowing that ip, 
/d, and do and aj obey bosonic commutation relations. We obtain the quasiparticle 


commutation relations: 

<P f (r ), <p f ( r') 


0(r),0 t (r') = <5(r - r') - C(r)C*(r')> 0( r ), 0( r O 


= 0, 


= 0. In the Bogoliubov limit the condensate wavefunction has no 

operator part, in which case ip may be written as ip(r) = T(r) + <p{ r). The condensate 
wavefunction has well defined phase and particle density and so may be expressed as: 
T(r) = a/ N 0 /A e* S( - r) p(r), where A is the area covered by the planar condensate. Note 

that the amplitude is normalized as A~ l f d 2 r p(r ) = 1. With these definitions, the usual 
bosonic commutation relations become: (p( r), (pfl r') = e lS ^ e ~ lS ( r ') <f(r, r'), where 

5(r, r') = <5(r — r') — A~ l a/ p(r) a/ p(r'). 

Next, we transform to the new Hamiltonian defined by K = H — pN = 
H — // Jd 2 rip' ip, then expand through second order in the operator part eliminating 
the linear terms by forcing the condensate wavefunction to satisfy the constraint 
(Ho — p, + g |T| 2 )'I' = 0. We arrive at the Bogoliubov Hamiltonian K = K 0 +K 2 , wherein 
zero-order and second-order operator terms are grouped into Ko and K- 2 respectively. 
These are defined as 


Ko= / drr T* 


Ho -h+ | I^W| 2 J 'h(r), 


K 2 = / d 2 rft( r) [H 0 - p + 2 g |T(r)| 2 ] </>(r) 


9_ 

"2 


d 2 r |[T*(r)] 2 0(r)0(r) + 0^(r)0'*'(r) [^(r)]" 


(31) 


Note that in addition to the kinetic operator we also have an arbitrary external 
potential in the first two terms, which in our case will be the periodic potential of 


the optical lattice. Equation (31) is quadratic in the field operators and so may be 


Uj{V) 




v j W 


d f 

3 


diagonalized with the appropriate field redefinition. To diagonalize Eq.(|3l|) we first 
apply the linear transformation <p( r) = e lS(r ' ) 

e -iS( u*(r) a) - vj( r) 


a, 


and <p\ r) = 

where the prime notation on the summation sign 
indicates that we are omitting the condensate from the sum. The afls and cc-’s inherit 
standard bosonic commutation relations from <p and (p\ and the spatially dependent 


transformation coefficients Uj(r] 


and Vj( r) obey the completeness relations 


^ [iij(r) u*(r') - »*(r) »,(r')] =i(r,r') 


( 32 ) 
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/ 


Y M r ) v Y) 

3 

o 

53 

* 

1 

(33) 

/ 

E K*( r ) 

- Vj(r) M*(r)] =0. 

(34) 


3 


So far, our discussion has taken place in two continuous spatial dimensions 
constrained only at the boundary by a trapping potential. We now want to translate 
to a formalism that fits a two-dimensional periodic optical lattice potential with 
honeycomb geometry. This is done by assuming a tight-binding limit at each lattice 
site. Formally, this corresponds to expanding the wavefunction in terms of a Wannier 
basis: functions which are localized and centered on each lattice site. The nearest- 
neighbor approximation then allows for a decomposition of the condensate and operator 
parts in terms of individual sublattices labeled A and B. In this new basis, the spatial 
dependence of the condensate and quasiparticle functions follows 


T(r) = ^ e* k (r rA ' > y/nA i e lSA i w(r — r^) + ^ e lk ' (r rB ^ y/n B . e lSs * w( r — r B ), 

A B 

/ 

0(r) = e lS(r) [ u jM r - r A)&j - v* j:A .(r - r A )a] 

A3 

/ 

+ e ' SM Y ( r - r n) Pi ~ v 1 b, (r — r B )yt . 

bj 


(35) 


(36) 


4-1. First method: Tight-binding limit followed by diagonalization of quasiparticle 
Hamiltonian 


We substitute Eqs. (35)-(36) into the Hamiltonian, Eq. (31), then take the long- 
wavelength limit while translating the exponential (crystal) momentum factors to 
coincide with the Dirac point. The continuum limit effectively converts the sublattice 
sums into integrals. By performing one of the integrations, over the A sublattice, say, 
while adhering to nearest neighbor overlaps, we obtain the affective Hamiltonian for the 
condensate and quasiparticles II = A 0 + K 2 where 


K n = / d 2 r 


ihcif:* A (r)V'if B (r) + ihciip* B (r)V* if A (r) + ^ |^(r)| 4 + ^ 


,(37) 


I<2 = Y^ / (j2r {wAi hci v iA r hci Vj.B Vvf A 

i.fc ^ 

+ 2U ajal v jA \if A \ 2 v* kA +2U Ml v jjB \fj B \ 2 v* kB 

~ 2 U \*Pa\ 2 aj&l (u j}A v* kA + u* kA v jA ) --U IV’bI 2 Ml ( Uj, B v* KB + u* k>B v jiB ) 
+ a]fl k hci u* A V*u k)B + (3ja k hc t u* B Vu k)A 
+ 2U a)a k u* A \^a\ 2 u k ^ A + 2 U fl]fl k u\ B \ip B \ 2 u k)B 

~\ U IVa | 2 ttjdfc {v*j,A U k,A + V kyA U* A ) - ^ U IV’bI 2 fljflk {v*j,B U k,B + V k}B U* j B ) 
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Sj bk hei vj 1 4 D a k B flj ex k hci Vj^B 'Tdii k ./\ 

- 2Ua j a k v jjA \^A\ 2 u kt A - 2 U fljfl k v jtB \^b \ 2 u k)B 

1 2 1 2 x x 

+ - |0a| OljaflUj^A Uk,A + V k! AVj,A) + 2 U Wb\ fljflk {Uj,BUk,B + V k} BV j}B ) 

- a)fllhc lU * A V*vl B - 

-2U a)oi\ u* A \ip A [ 2 v* kA - 2 U fl)fl\ u* B |0.b| 2 v* kB 

+ \u |0a| 2 oi]al ( u* }A u* kA + vl A v* A ) + ^U \^ B \ 2 « B < B + v I,b^b) 


■ (38) 


Here we have defined the condensate two-spinor in terms of the A and B sublattice 
components T(r) = [-0 a( r ), 0s(r)] T , and the Dirac operator is defined as V = d x +id y . 
Next, we isolate the first six terms (terms with the daggered operator to the right) and 
write them as a matrix contraction of two pure operator valued vectors 




Da,b \ f d \ 


(39) 

(40) 


where 


A u , v = 2 XJv jtA |0a| 2 v* kA - -U |0a| 2 (u jA v* kiA + K,a v j,a) > 
Bu,v = ZUv j>B |0b| 2 V* k B - \u |0b| j (u j>B v* k)B + u* k B v^B ) , 


D a ,b = hciVj : A'D*v k B , 
D Bi a = hciVj^V v* k A ■ 

The eigenvalues are then obtained by 
A lf tj A D a ,b 

>B,A 


det I 


lU,V 

D f 


B h 


= 0 


(A UjV — A) ( B UjV — A) — D a ,b D b ,a — 0, 

A± = ^ ~ ^ 2 \!(^“^ — Bu,v) + 4 Da,b D Bi a , 

and the corresponding eigenvectors follow 


H + = 


1 

Db,a 

(-^rt Bu,v ) 


The unitary matrix that diagonalizes Eq.(39) is 


u = -p 

V2 


1 

Db,a 


1 

Db,a 


D+ B u , v j (A_ Bu, v ) 


The first six terms ill Eq. (38) may be expressed ill the new basis as 
A+{jt} c+y c + f. + A_j,r} 0-,j l~ t k i 


(41) 

(42) 

(43) 

(44) 


(45) 

(46) 

(47) 


(48) 


(49) 

(50) 
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where we have included the j, k subscripts on the eigenvalues to be fully descriptive. 
The new quasiparticle operators can be written in terms of the old ones as 

D* 


r f - 
C ±J ~ 


1 

71 


a] + 


B,A 


t 

B* ) r ' j 

u.v) 


p 


(51) 


Note that the right hand side is /c-dependent which is implied on the left. The substance 
of the transformation is contained in the momentum and space-dependent eigenvalues 

(52) 


Km = Uv j)A \^a\ 2 v* kA - \u |^ A | 2 ( u jtA v* KA + u* kA v jjA ) 


+Uv j)B \^b\ 2 V* k B - - U I^bI 1 {u j)B V* KB + K,B V j,B ) 


± 


Uv j)A \^ A \ 2 v* kA -\u |^ a | _> (■ u jiA v* kA + u* kjA v jjA ) 


UVj,B I'lpBl 2 V* ktB + \u IV’bI 2 {Uj, B V* KB + U* kjB Vj,B) 


I 2 


+ (hcif v jA (V*v^ B )v jiB iVvl A )} 


1/2 


The next step is to constrain the quasiparticle amplitudes in Eq. (53) (the u 1 s and 
v’s) in order to diagonalize the Hamiltonian with respect to the momentum indices j 
and k. First, we let 

hciv jA V*v* kB = 2Uv jA IV’aI 2 ^ - \ U |^a| 2 (■ u jiA v* kA + u* kA v jA ) 
hciv jA Vv* kA = 2 Uv j>B |Vib| 2 vIb ~ \u IV’sl 2 {u jA v* kB + K,b v 3 ,b) , (53) 


and then substitute these into Eq. (53), which reduces the two eigenvalues to 

Km = ~ h v j, A v t, A + 2 Uv jA \^ A \ 2 v* kA ~\u |V^a| 2 (■ u jjA v* k}A + u* kA v jA ) 

~hVj,BV* k)B + 2 Uv jjB \^b\ 2 V* ki B --zU \^b\ 2 { u j,B v k B + K^V^b) , 


and 


X -{jk} 0, 


(54) 


where we have reinserted the chemical potential terms. It is important that Eq. (53) 
depend only on one index so that quasiparticle amplitudes for different eigeneneregies 
are not coupled. Dividing Eq. (53) through by Vj, A and Vj iB , respectively, cancels all 
j-index terms except for ones that appear as Uj A /v ]A and u j B /v J)B . To completely 
decouple the j-k modes, we must ensure that Uj A /vj A = Uj A /vj )B = h( r )> be., the 
amplitudes for any given quasiparticle mode have the same relative spatial form. We 
can then rewrite A +{jk} as 

Km = (55) 

1 1 o 1 2 

2 hc i V 3 ,A V * V k,B - 2 A* v i,A v*k,A + u V J,A \ipA I v* KA --U\if> A \ (u jA v* kA + U* kA V jA ) 

1 1 2 1 2 

+ 9 hc l V 3,B V V* k A - -ll V jB V* kiB + UVj t B W’B \ V* k B - - U \^ B \ {u jA V* k B + U* k B V jA ) . 
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Finally, we impose the constraints 

- 7 E kV*k,A = 7 hc i v * v k,B - tKa + b* 7 l^l 2 _ 7^ l^l 2 


^ ^ ^ 6 «,-D 2 Ir 1 ^ 

1 1 1 1 2 1 2 

-EjV jA = -hc{D‘Vj,B - + ^U\il A \ Vj,A - j U\i> A \ Uj,A- 


(56) 

(57) 


Multiplying Eqs. (56)-(57) by v hA and v kA , respectively, and using the property that 
f d 2 r Vj ,b T> v k A — J d 2 r{V*Vj ) B ) v* k a > we ma y separate out 1/2 of each derivative term 
in Eq. ((561), which reduces the non-derivative terms in the first, line of Eq. (56) to 


— + Ej)vZ >A Vj iA ■ 


(58) 


We may reduce the second line using the other half of each derivative term, thereby 
condensing the eigenvalues down to 


A + {« = ~A E * + E j)( V k,A V J.A + B Vj.B ) ■ 


(59) 

The next six terms in Eq.([38j) may be diagonalized in a similar way yielding the 
eigenvalues 

7 +{i^l // ^k,A T 2 U aj A ^ k ,A 

1 2 

--U\ij) A \ ( V* jA Uk,A + v k)A u* A ) 


-!XU* B U ki B + 2 Uu* jiB \i’ B \ 2 u kt B 

1 2 

--U^b] (v* B u k B + v k)B u* jB ), 


and 


^-{jk} — 0 • 


Following our previous steps, we obtain 


^+{jk} . (E k + Ej ) (u k A Uj,A + u k: B U j,B ) • 


(60) 

(61) 

(62) 


Combining Eqs. (59) and (62), and inserting the quasiparticle operators, reduces the 
first twelve terms in Eq. (38) to the expression 

1 \ ' f T 1 

2^2 d 2 r (Ej + Ek) c ] + jC +ik (u* KA Uj,A + u* kB u jtB ) ~ c+^c ] +k (v* kA v jA + v* kiB v jiB ) • (63) 


j,k 


For the special case where j = k, we may further combine the terms at the cost of an 
extra c-number term to arrive at 
1 x ( f 

~ g / d2 ‘ r 2E k( V *k,A v k,A + V* k B V k , B ) 

k J 

+ 7 / d2r ( E i + Ek ) d +,j d +,k( U k,A U j,A - V*k, A Vj,A + U* k>B U jiB ~ V^ B V j)B )- 

A U J 


( 64 ) 
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Applying the completeness relations Jd 2 r {u* k A Uj,A ~ a v PA)—^'c'j an d / d 2 r ( u* k B Uj tB ~ 
v t A v j,s)~^i,v contracts Eq. (64) down to 


E f ^ rE k ( \ v k,A | 2 + \Vk,B?) + E E ‘ 

1. J 1- 


at 


k C +,k C +,k ■ 


(65) 


Diagonalizing the rest of Eq. (38) (terms with no daggered operators and ones with only 
daggered operators) by capitalizing on the j-k symmetry of terms such as J d 2 r Uk,AVj,Ai 
and anti-symmetry of the (.Ej — Efl) factor, we obtain the final form of the interacting 
Hamiltonian 


H= jd 2 r [ihci if}* A (r)(d x + idy)^ B ij) + ihci^* B (r)(d x - id y )i/j A (r) (66) 

+ f//2|^(r)| 4 + D/2|V^(r)| 4 ] 

’ r ' 

- E E > I d2r v J* v i + E E i d +’i > ( 67 ) 

3 3 

with the resulting constraints on quasiparticle amplitudes given by 

hciV*u jjB + (2 U |-0a| J - At) Uj iA - u |-0a| J v jt A = E jUj:A , (68) 

hciVu jjA + (2 U IV’sI 2 - d) u j,B ~ U l^sl 2 v j>B = EjU jtB , (69) 

-hciVv j:B - (2U \^a\ 2 - d) v jt A + u\fl> A \ 2 u j: A = E jVjiA , (70) 

-hciV*v jA ~ (2 U |-0s| 2 - d) v j,B + u \ip B \ 2 u jt B = EjVj }B . (71) 


4-2. Second method: Diagonalize quasiparticle Hamiltonian then impose tight-binding 


Although the first method is cumbersome it is the more rigorous approach and instils 
confidence in the final constraint equations. A shorter approach is to first obtain the 
usual Bogoliubov equations for a condensate not confined in a lattice, and then apply 
the tight-binding limit directly. The Bogoliubov Hamiltonian is 


H= / d 2 r T* 


Ho -d + | I^WI 2 ^(r) 


E S i d\\ v flr)\ 2 + J2 E . 


jOljOlj , 


with the constraint equations (BdGE) given by 


Cuj — g |T| 
C* Vj - g 


Vj EjUj 


l *| 2 “j = -EjVj 


(72) 


(73) 

(74) 


In Eqs. (73)-(74), C is a differential operator containing terms that couple the 

An additional implicit constraint is that T 
To pass to the tight-binding limit we 


quasiparticle and condensate velocities, 
satisfies the nonlinear Schrodinger equation. 


express all spatial functions in Eqs. (72)-(74) in terms of Wannier functions for the 


individual sublattices, and evaluate the Bloch plane wave factors at the Dirac point 
momemtum. Adhering to nearest-neighbor overlap for on-site Wannier functions, we 
integrate out spatial degrees of freedom (which splits the honeycomb lattice into A 
and B sublattices), regroup terms into finite differences, and then take the continuum 
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limit. Equation (72) then transforms to Eq. (67), while Eqs. (73)-(74) transform to 
Eqs. (|68|-([7l| with several additional derivative terms contained in C as follows 


- [V 2 + iV 2 0 + 2 iVf V- (V0) 2 1 Uj , (75) 

2m 

where 0 is the condensate phase. After going through the steps that culminate in the 
tight-binding continuum limit, these terms transform to 


ih Cl V*u kM A) + [-iftTiV0A(B) ■ V + &r 2 |V0A(B)i - ihr 3 (V 2 0A(B))] u k)A (B ), (76) 

where the coefficients encapsulate the spatial integrals as follows: Tioc | [dr 
t 2 oc | fdr w* a {S 7 (fl)w a\-, T 3 oc JdrV0| 2 w^- These extra terms depend on the 
condensate phase 4>a(b), and so couple the superfluid velocity to the quasiparticle 
excitations. In particular, the term with coefficient T\ depends on the direction of 
quasiparticle emission relative to the motion of the condensate. The relativistic linear 
stability equations, Eqs. @)-0), may be expressed in compact notation as 



where 

A = U diag(|V>n| 2 , |0n| 2 ), (78) 

Ak = A k ■ 1 2 , (79) 

[D] 1A = -/jl + 2U 10a! 2 - ihrflJ^A ■ V + hr 2 |V0a| - ihr 3 (V 2 0^) , (80) 

[D] 22 = ~n + 2U - ihr{V(f>B ■ V + hr 2 |V0 S | - ihr 3 (V 2 0 B ) , (81) 

P]i,2 = [^2,i = iteiV', (82) 

and 1 2 is the 2x2 unit matrix. 


5. Stability of vortex solutions 


Two independent derivations of the RLSE in Sec. [4] and their reduction to the BdGE, 
which we discuss in Sec. [6j establishes Eqs. (77)-(82) as the correct method for computing 
the low-energy structure (quasiparticle states and eigenenergies) for arbitrary vortex 
solutions of the NLDE m Radial profiles for vortex solutions of the NLDE are plotted 
in Fig. [3j with details of our solution methods presented in [43], Briefly, the plots in 
Fig.0 were obtained by solving the NLDE in plane-polar coordinates for the case of 
cylindrical symmetry, which read 


- hc t (d r + f B (r) + U | f A (r)\ 2 f A (r) = /x/ A (r), (83) 

hc t (d r + + U \f B (r)\ 2 f B (r) = fif B (r) . (84) 

for the upper and lower two-spinor component radial functions f a and f B , and £ G Z the 
angular quantum number coming from imposing single-valuedness of the wavefunction. 
The most immediate and pragmatic concern is the combined effect of the honeycomb 
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Figure 3. (color online) NLDE vortex radial solutions, (a) Bessel solution for £ = 3; 
(b) numerical solution for £ = 4; (c) vortex/soliton; (d) Anderson-Toulouse vortex; (e) 
ring-vortex solution for £ = 4; (f) ring-vortex/soliton solution; (g) Mermin-Ho vortex; 
(h) half-quantum vortex. In each plot, the upper and lower spinor radial solutions are 
indicated in red and blue, respectively. 


lattice geometry and the inter-particle interaction on the lifetime of a vortex. It should 
be emphasized that the presence of an infinite tower of negative energy states below 
the Dirac point seems to imply that a condensate residing there will eventually decay 
provided there is a mechanism for energy dissipation into noncondensate modes (i.e., 
secondary interactions with thermal atoms) J Generically, negative energy states are 
present for moving condensates for which excitations subtended by a backward cone have 
negative frequencies [SO] • Moreover, when a vortex is present small displacements of the 
core from the symmetry axis of the trap results in a precession of the core, which, when 
combined with dissipation, causes the vortex to spiral to the edge of the condensate. 
In the absence of a periodic lattice potential this dynamical process is known to be 
driven by the anomalous modes in the linear spectrum, i.e., modes with negative energy 
and positive norm [83] also called Goldstone modes. The time for a vortex to spiral to 
the edge of the trap would then define its lifetime. In the absence of the lattice this 
precessional motion is canceled by introducing rotation to the trap [83] [ST] , a result 
which we suspect to be true in the lattice case as well. 

To undertake a full treatment of the lifetime would mean computing this spiraling 
time and then comparing it with the lifetime that we compute here due to the dynamical 
instability from the complex frequencies. The lifetime of the vortex would then be the 
smaller of the two values. Nevertheless, in cases where dissipation is weak and the 
vortex is centered on the symmetry axis of the trap, the dominant source of instability 
arises from the complex eigenvalues associated with RLSE modes. We will limit our 
analysis to the effect of the latter, and regard the negative real part of the eigenvalues 
from a standpoint of metastability. Physically, the complex eigenvalue gives rise to 

f We remind the reader that this infinite tower of negative energy states is only in the Dirac cone 
approximation. 
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fluctuations in the angular rotation of the vortex spinor components [ZB. In the case 
of the NLDE this is a result of internal “friction”, i.e., energy exchange, between the 
two spinor components displayed in the complex derivative terms of the Dirac kinetic 
energy. This drag force between the two vortex components (or between vortex and 
soliton) eventually causes substantial degradation of the vortex. This is the measure 
that we will use to compute vortex lifetimes. 


5.1. Numerical solution of the relativistic linear stability equations and vortex lifetimes 
The stability of a particular condensate density and phase profile such as a vortex 


is arrived at by expanding Eq. (77) and expressing differential operators in terms of 


suitable coordinates, for example polar coordinates for a vortex, then using separation 
of variables for the quasiparticle amplitudes with the appropriate form of V’A(B), be., 
solutions of Eqs. (83)-(84). This yields a set of first-order coupled ODE’s in the 


radial coordinate to be solved consistently for the functions Ua(b){ t )i Va(b )( t ) and the 
eigenvalues E We discretize the derivatives and functions using a forward-backward 
average finite-difference scheme, then solve the resulting discrete matrix eigenvalue 
problem using MATLAB function Eig. In Fig. [4] we have plotted the real and imaginary 
parts of the first 20 eigenvalues, labeled by the quantized quasiparticle rotation number 
n G Z, for the vortex/soliton solution which we discuss in previous work |J5]. The 
lowest modes are anomalous with negative real parts and positive, nonzero but small, 
imaginary parts. 
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Figure 4. (color online) Anomalous mode frequencies for the vortex/soliton. The 
real part of the anomalous mode frequencies are plotted in (a), the Imaginary parts 
are plotted in (b). The horizontal axis labels the excitation mode determined by the 
quasiparticle angular momentum quantum number n. 

Convergence of RLSE eigenvalues for the l = 1 vortex/soliton background as a 
function of the grid size N used in the 4 A" x AN matrix problem is displayed in Fig. [5j 
where we have plotted the real and imaginary parts of the eigenvalue for the lowest 
excitation mode. The lifetime of a particular vortex solution can be computed by 
examining the lowest quasiparticle rotation mode n — 1 , since at very low temperatures 
this mode dominates the spectrum. The lifetime is then characterized by the reciprocal 
of the imaginary part of the associated eigenvalue, i.e., lifetime = H/lm ( E_i ). Here, the 
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—1 subscript refers to quasiparticle rotation relative to the vortex rotation. Eigenvalues 
for the lowest quasiparticle rotational mode and the associated lifetimes for all of onr 
solutions are listed in Table [U 



Figure 5. (color online) Convergence of RLSE for the vortex/soliton. Real (a) 
and imaginary (b) parts of the lowest anomalous mode. The horizontal axis shows 
the number of steps and the energy of the lowest excitation of the vortex/soliton 
corresponding to n = 1 in Fig. [4] is plotted on the vertical axis. 


Solution type 

Quasiparticle energy [nK] (±10 4 nK) 

Lifetime [s] 

Complex topological vortex 

2.231 - 4.174 x 10 2 * 

1.969 x 10“ 5 (±4 x 10” 12 ) 

Topological vortex 

8.184 x 10' 3 - 1.066 x 10 3 i 

1.941 x 10“ 5 (±6 x 10" 13 ) 

Ring-vortex 

-4.181 - 1.599 x 10~ 2 i 

0.5295 (±3 x 10~ 3 ) 

Ring-vortex/soliton 

-4.203 + 2.022 x 10' 3 * 

4.043 (±2 x 10- 1 ) 

Vortex/soliton 

-4.211 ±2.141 x 10' 3 * 

3.841 (±2 x 10- 1 ) 

Mermin-Ho vortex 

2.818 x 10 2 ± 1.066 x 10 5 i 

7.712 x 10^ 8 (±7 x 10- 17 ) 

Anderson-Toulouse vortex 

-4.202 ± 2.033 x 10" 3 i 

4.041 (±2 x lO^ 1 ) 

Half-quantum vortex 

2.818 x 10 2 ± 1.066 x 10 5 i 

7.712 x 10" 8 (±7 x 10- 17 ) 


Table 2. Stability properties of NLDE vortices. Lifetimes are computed using the 
value of the interaction U in Table [l] and the formula lifetime = h/lm. {E_\). The 
various vortex solutions of the NLDE are plotted in Fig. [3] and derived in detail in |43j . 
Note that solutions with similar boundary conditions have lifetimes of the same order 
of magnitude. 


To understand the character of the quasiparticle modes we must consider the spatial 
functions associated with each eigenvalue. Radial profiles for the n — 1 quasiparticle 
excitation in the vortex/soliton background are shown in Fig. [b] They are bound states 
near the core of the vortex localized specifically at the point where the soliton and vortex 
components of the background are equal (see Fig. |3](c)). Physically, the imaginary part 
of the eigenvalues imply a transfer of energy between the vortex and soliton components 
through quantum fluctuations. In particular, each component acquires quantum 
admixtures from different rotational modes as well as local shifts in amplitude from phase 
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Figure 6. (color online) Plots of lowest quasiparticle excitation for the vortex/soliton 
configuration. Used with permission inl¬ 


and density fluctuations, respectively. Mathematically, the full quasiparticle operator 

r „ „ -| T 

with time and spatial dependence for this mode is 4>{r,t) ~ <f) B _i(r,t ) , 

where the quasiparticle spinor operators are 


4> A ,-i(r,t) = e ^ t/h e t6 u A _i(r) d_i + e zE ~ lt/h e l0 v* A _flr) , (85) 

fe-i(r,t) = e~ iE ~^ h u B ^ (r) + ^-^^(r) • (86) 


As discussed previously, relative to the vortex background the quasiparticle has rotation 
£ = —n = — 1, which has the effect of reducing the rotation of the vortex. Note that the 
expression for the operator 0(r, t ) is approximate since we have truncated the sum over 
quasiparticle modes after the lowest mode. We recall that the spatial functions have the 
properties u A -i(r), mb,-iW ~ 1CU 2 and v A -flr ), v B _i(r) ~ 10 -5 (see Ref. mi). where 
all are peaked in the “notch” region ^Dirac < r < 2^Dirac, and where the absolute values 
of the slopes of the soliton and vortex are maximum. In this region, the normalization 
integrals (one for each sublattice) are given by 

Jd 2 r [\u A{B ) i(r)| 2 - |n j4(B )_i(r)| 2 ] > 0. (87) 

This combination of positive norm and negative Re(E_!) signals the presence of the 
anomalous mode. In Sec. [6j we will see that these bound quasiparticle modes solve 
the Majorana equation, which predicts an additional zero energy mode localized at the 
same distance from the center of the vortex. 


6. Connection to other theories 

In this section we examine several reductions of the RLSE to other equations familiar to 
BECs, superconductivity, graphene, and high energy physics. Our results demonstrate 
the variety of substructures contained within the RLSE framework. Note that we adhere 
to the weakly interacting regime through all of our derivations as explained in Sec. [2j 
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Mappings of the RLSE to other more fundamental equations proceeds by resolving the 
RLSE solution space into a parameterization with respect to two measures: 1) the ratio 
of the on-site energy, which we denote t 0 , to the chemical potential p\ and 2) the ratio 
of the quasiparticle kinetic energy, cip, to the boson interaction strength U. The on-site 
energy is calculated as the average value of the Laplacian plus lattice potential over 
same-site Wannier functions. As such, t 0 encodes the sublattice energy offset into a 
uniform mass gap at the Dirac points of the continuum theory. Changing to/p tunes 
the spectrum between a linear (relativistic) and quadratic (nonrelativistic) dispersion. 
In contrast, changing Cip/U from large to small values tunes the spectrum from a pure 
particle-like dispersion to one characterized by an equal admixture of particles and holes. 
The dimensionality of solutions to Eq. (77) experiences a corresponding change over the 
parameter space from a single-component Schrodinger-like solution (for to/p, cip/U ~ 1) 
to a four-component spinor solution (for to/p, Cip/U <C 1), the latter similar to Nambu- 
Gorkov states in superconductors arising from doubling of Fermion degrees of freedom. 
Thus, the quasiparticle spectrum is spanned by a two-dimensional parameter space 
highlighting the similarity between Bogoliubov and relativistic structures. 

To quantify our discussion we look for a non-relativistic reduction of Eq. (77) 
by working first from the lattice form of the NLDE since the hopping terms are the 
same as those for the RLSE. We recall that the standard massive Dirac equation has a 
well defined non-relativistic limit to the Schrodinger equation [MJ. The proof uses the 
fact that in the low-energy limit the mass term, proportional to me 2 , is the largest 
contribution to the energy. In particular, the two-spinor formulation of the Dirac 
equation is comprised of two coupled equations. The procedure involves isolating the 
mass term in one equation and substituting the resulting expression into the second 
equation. The substitution converts the first-order spatial gradient to a second-order 
Schrodinger kinetic term for small relative kinetic to mass energy. This effectively pushes 
the mass dependence out into smaller correction terms which may be neglected. Similar 
steps may be implemented in our case but we must first introduce an offset between the 
sublattice potential well depths (a mass gap) so that we obtain the desired curvature in 
the spectrum near the Dirac points, effectively introducing a non-relativistic regime. 

Starting from the discrete NLDE for a single Dirac point and following similar steps 
as in our previous work [27], we obtain 


pfl>A j = ~th {i’Bfl 


zk-^3 


+ ^B j . ni e 


zk-5i 




) - t 0 if Aj + u | if Aj 


12 


( 88 ) 

(89) 


+ 'Pb j . 

= -th (^A s e~' kl ‘ + *pA J+ni e~ tkSl + Vt i+ „ 2 e~* k "’ 2 ) + t«VB, + U\tp Bl \‘ipB J 
where th, to, U, and k are the hopping, same-site, and on-site interaction energies and 
crystal momentum, respectively. The <5’s, n’s, and 2D vector indices j indicate the 
lattice vectors described in Fig. [7] In Eqs. (88)-(89), t 0 is the sublattice offset equivalent 
to a spectral gap 2|£o|- For weak interactions, the on-site energy can be made much 
larger than the contact interaction strength by tuning the lattice potential so that 
|p ± t 0 1 >> U. After inserting the correct values for the lattice vectors and solving 
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B 



c) 



Figure 7. (color online) Characterization of a honeycomb lattice, a) Hexagonal lattice 
structure, b) Nearest neighbor displacement vectors, c) Reciprocal lattice. Used with 
permission m- 


Eq. (89) for ip Bj , to zeroth-order in U/\fi — f 0 | we obtain 

t h 


fpB, = 


h ~ to 


bPA t +^/ 2 ' /3 +^,e-‘ 2 ’' 3 ) 


(90) 


From Eq. (90) we may write analogous expressions for neighboring sites by shifting the 
indices using the lattice vectors n,-, 


lp B . = 

T Uj— ni 


th 


h — to 
th 


+ 'ipA-e 12 '*/ 3 + ip a-,, ,< 

— Til ' ^3 ' + n l) 


'll) b • — 

u J=>j-n 2 _ -L 

fl TO 


-j-n 2 + 1pA,_ 


J2ir/3 




+ fl’A.e' 


i27r/3j 

(91) 

i27r/3\ 

(92) 


Substituting Eqs. (|90|)-(j92|) into Eq. (|88|), expanding complex factors and regrouping 
tc 

tl 


the terms to form finite differences, we arrive at the expression 
t 2 

l h 


= 


{(dPj+n i 2 ^j + 'lpj — ni ) + (ipj-\-n 2 2 Ipj -\- lpj— n 2 ) 


J 2(/i —1 0 ) 

"h ( 1pj+(ri2—n\) 2 Ipj + Vh'—(n 2 —m) ^ iV 3 [(Vh+ni Ipj) + (Ipj lpj— ni ) (fPj+ri2 V(?) 

{rfj ~ V’i—raa) “I" (fPj+(ri 2 -ni) ~~ V’i) + (V’j — V’j—(n 2 —ru))] } — t$lpj + £/ |^j| Ipj ■ (93) 


Equation ( |93| is a discrete nonlinear Schrodinger equation for the honeycomb lattice in 
the sense that it has as its continuum limit the usual nonlinear Schrodinger equation with 
cubic nonlinearity. Substituting the correct continuum forms for the finite differences 
and then expressing the result in rectangular coordinates, we obtain 


hi’ = 

t 2 a 2 \ (3 d 2 1 8 2 

2(fi — t 0 ) + Ady 2 


d 2 ip ,\/3 / r-dip dip 
dy 2 a \ dx dy 


VS d 2 \ f3d 2 - 

2 dxdy J dx 2 

/-dip dip n dip\ 
^Sx ~ ~dy + 2 dy) 


1 d 2 V3 d 2 \ 
4 dy 2 2 dxdy J 

- t 0 ip + U \ip\~ ip, 


ip 


( 94 ) 
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which finally reduces to 


(/i + t 0 )ip = - Clh , V 2 il> + U \i/j\ 2 -if ), 


(95) 


(h - to) 

where we have substituted in the effective speed of light q = \f+t h a/2fi. Performing the 


same steps with Eq. (89) gives a second equation 

c? 


(h - to) = 


2 h 2 


(/I + to) 




(96) 


Next we examine two limits. For t 0 <C p, Eqs. (95)-(96) describe two propagating 
modes where the effective mass and total energy are of the same order. To lowest order 
in t 0 /p Eqs. (95)-(96) become 

(97) 

(98) 


and 


p 2 ^ + c 2 h 2 V 2 ^ + pt 0 fl> — pU \ip\ 2 ip = 0 , 
p 2 ^ + c 2 h 2 V 2 ^ — pto^ — pU |^| 2, 0 = 0. 


Reintroducing time dependence through p 2 —> —K 2 d 2 and dividing through by c 2 h 2 gives 
1 d 2 ^ „ 0 , m 2 c 2 , U 2 


:? dt 2 


1 d 2 fl> 
c? dt 2 




me 

—fl 2 


V’ + -2^IV’l'V’ = o, 


m 2 cj 


__,9 . m u , 

V 2 V> + + 


c 2 h 2 

u 2 

c 2 h 2 


V = o, 


(99) 

( 100 ) 


where we have defined the mass m = \/\pto\/c 2 and interaction strength U = \J\pU\. 
Equations (99)-(100) are nonlinear Klein-Gordon equations describing a tachyon mode 
with imaginary mass in Eq. (99), and an ordinary Klein-Gordon mode with real mass 
in Eq. (100). In contrast, if we tune the lattice potential offset so that t 0 ~ p the 
mode described by Eq. (95) has a very small effective mass and large energy, whereas 
the mode in Eq. (96) will have a very large mass and small energy. In this case the 
mode in Eq. (96) gets “frozen out” and we are left with only one propagating mode in 
Eq. (95). Here multiplication by the total energy p+to does not cancel the effective mass 
p — t 0 in the denominator of the gradient term. Reintroducing the time dependence by 
p + to —>• ihdt and the effective mass m = (p — t 0 )/2c 2 , Eq. (95) reduces to the nonlinear 
Schrodinger equationt 


ih^ +-^-V 2 ^ - U \^\ 2 ^ = 0 . (101) 

Thus, tuning to interpolates between a Dirac and a Schrodinger structure with Klein- 
Gordon bridging the two. One may understand the intermediate Klein-Gordon result 
through a general argument by noting that any reduction of the Dirac equation to the 
Schrodinger equation must modify both the single-particle dispersion as well as the 
spin attached to each excitation mode. To be precise, two regimes are identified: one 
associated with binding two spin-1/2 modes into a single spin-0 mode at lower energy 


f This step can be justified formally from the Heisenberg equation of motion for the wavefunction 
starting from the operator formalism, but such justification is well known from theory of NLSE. 
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resolutions (Eqs. (88)-(89) to Eqs. (99)-(100)), and one associated with a crossover 


from relativistic to classical dispersion (Eqs. (99)-(100) to Eq. (101)). This observation 


applies not only to the fundamental case, but also in the honeycomb lattice picture with 
regards to pseudospin. 


Applying the same steps to the RLSE, Eq. (77) yields the BdGE 

— A r , 


u) 

1-/ 

l • i 

1 \ 


-^V 2 + A r 


2m 


A 


ALy 2 

2m v 


A, 



( 102 ) 


with A m = —p + 27/|0| 2 , A p = U |0| 2 , where 0 is the condensate wavefunction for 


either of the decoupled sublattices and the effective mass is the same as in Eq. (101), 
m — (p — to)/2c 2 . Note that we have suppressed explicit space-time dependence 
in Eq. (102) for clarity. In the particle regime for large characteristic momentum, 
q|p| U, the particle and hole amplitudes satisfy u v and Eq. (102) reduces to 
the standard Schrodinger equation for a particle moving in the potential V = — A m . 

Next, we look at the case of single-mode approximation for the pseudospin degrees 
of freedom in Eq. (77), i.e., where the sublattice backgrounds are equal 0^ = 0# which 
also implies that ua = ub = u(r,t) and va = vb = v(r,t). One then finds that the 


system Eq. (77) reduces to the Andreev equation 



u) 

1 -/ 

V v ) 

1 \ 


-ihcip ■ V + A, 
A„ 


-A ? 

ihcify ■ V 



(103) 


The unit vector p in Eq. (103) points in the direction of quasiparticle propagation. 


Here we have chosen the case of zero background flow \74>a,b = 0 as we will do for 
the remainder of this section except for the vortex background. Equation (103) is 


the Andreev equation for propagation through a medium comprised of both normal 
and superconducting regions [55] . The spatially dependent pairing and mass terms 
are A p (r) = f/|0(r)| 2 and A m (r) = 2T/|0(r)| 2 — p. In this analogy the condensate 
wavefunction 0(r) stands in for the order parameter in a superconducting medium. 


Equation (103) describes slowly varying particle and hole functions u( r) and v(r) 
split off from an overall rapidly oscillating plane wave portion which moves in the 
direction p. Thus, we should expect similar exotic scattering such as specular and 
retro-reflection IHoj . 

Next, we look at the particle regime where the particle component is dominant, 


ua(b) v a,(B)- bi this regime Eq. (77) reduces to the Dirac equation 

id y ) 



A A 

-ihci(d x + idy ) 


-ihci(d x - 
A b 



(104) 


where a potential term appears A^( S )(r) = 2U\'i/jA(B)( r )\ 2 ~ h- Equation (104) 


further reduces to the massless Dirac equation in the case of a constant background 

|^,B)| 2 =M/(2 U). 

Interestingly, zero-mode solutions (E = 0) of the RLSE occur as well and we 
find that these solve the Majorana equation which is implicit in the RLSE for certain 
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background configurations. To see this we set Ek = 0 in Eq. (77), which decouples the 


system into two sets of equations in the extreme long-wavelength regime characterized 


by \ua(b)\ = \va(b)\- In this regime Eq. (77) gives two copies of the form 



= 0 


(105) 


—ihcfld x — idy) 

-ihcfld x + id y ) A b 

where the potential terms are A a{b) = U\^a{B)\ 2 — h- For a uniform condensate, i.e., 
far from any vortex cores, the asymptotic choices are U\i/ja(b) | 2 —> /x, 0. In both cases 
Eq. (105) offers no solution. However, for the vortex/soliton there is a “notch” in the 


order parameter near the core, where \^a\ 2 = \^b\ 2 < h/U =>- A j4 ( S ) < 0, in which case 


Eq. (105) reduces to the Majorana equation 
— ia ■ Vi/j c + mi/) = 0 , 


(106) 


where 0 C = it])* = [u* A , u * B ], m = |A^( S )|, and 0 : [R 2 —> [R 2 . Equation (106) supports 


real solutions with linear dispersion and has been studied extensively in its original 
mathematical form EH and more recently in condensed matter physics intimately 
associated with topological insulators [55]. In their present incarnation these Majorana 
zero modes also occur in the core of nonlinear Dirac vortices with higher winding (t > 1 
in Ref. [35]) where both spinor components vanish |0 a(B)(O, 0)\ 2 = 0. In this case the 
mass term in Eq. (106) reduces to the condensate chemical potential m = ji. In the 


superfluid context the meaning of the Majorana zero mode is of a zero-energy pure 
spatial density fluctuation associated with rigid translations of the vortex core. Here 
phase fluctuations only appear as finite-energy fluctuations in the vortex rotational and 
translational motion. For the vortex/soliton the zero mode is a circular ring reflecting 
the symmetry under both rigid rotations as well as translations of the vortex. In Fig. [8] 
we summarize the various types of reductions of the RLSE indicating the conditions or 
limits for each equation type. 


6.2. Mapping to relativistic Bardeen- Cooper-Schrieffer theory 

In this section we discuss the modifications needed to connect the RLSE to relativistic 
BCS theory. Here we capitalize on an important property of the NLDE and RLSE. 
This is that repulsive interactions for bosons in the honeycomb lattice break the valley 
particle-hole exchange symmetry at the Dirac point in a significant way such that 
an additional sign change of the interaction restores the symmetry. More properly 
stated, the noninteracting theory is invariant independently under charge conjugation 
(C), parity inversion (' V ), and time reversal (T). Repulsive interactions break T and 
C, but the symmetry-breaking cancels in such a way as to preserve the full CPT 
symmetry [27]. Consequently, a parity inverted positive energy solution (valley particle) 
can be interpreted as a negative energy solution (valley hole) in a theory with attractive 
interactions but without parity inversion. Stated differently, a theory of particles with 
repulsive interactions is equivalent to a theory of holes with attractive interactions. 
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Figure 8. (color online) Schematic of reductions of the RLSE. Limiting theories for 
the RLSE are displayed for relative strengths of the lattice potential offset and chemical 
potential along the vertical direction and quasiparticle momentum versus interaction in 
the horizontal direction. The dashed lines indicate crossovers in the underlying lattice 
theory. 


To complete the mapping to BCS theory we introduce a mass term and nearest- 
neighbor interactions at the lattice scale to couple the different spinor components. The 
mass term is obtained through an asymmetry in the honeycomb sublattice potential 
depths, an intermediate step in populating Dirac points, as we have explained in [43]. 
The various types of relativistically invariant interactions may be constructed using 
nearest-neighbor interactions as follows. Specifically, the symmetry of the nonlinearity 
in the NLDE determines the symmetry of the superconducting order parameter and pair 
potential in the corresponding BCS analog equations E 3 - The vector-vector interaction 
can be obtained by including repulsive nearest-neighbor interactions. A scalar-scalar 
type coupling can be realized similarly, but by using attractive (instead of repulsive) 
nearest-neighbor interactions in addition to the repulsive on-site interactions. The 
spin and pseudo-spin symmetric terms are characterized by an alternating sign for the 
coupling between the two spinor components. This type of coupling may be realized in a 
lattice setting via Feshbach resonances using a beam with the proper spatial modulation 
to produce interactions whose sign alternates between neighboring lattice sites. Pseudo¬ 
scalar forms can be realized by eliminating on-site interactions while retaining repulsive 
nearest-neighbor interactions. 

The case of scalar-scalar coupling in the NLDE with equal on-site and nearest- 
neighbor interactions U = U nn and mass term m s tf (see ref. [43]) elevate the RLSE to 
the form 


u k 

Vk 


( 107 ) 













Superfluid fluctuations and emergent theories 


35 



ihcia ■ V + m s c 2 • 1 2 + gcr M A M —iA p a y 

iA p a y — ihcia ■ V + m s c 2 • t 2 + 

where A p (r) = f/[|'0 J 4(r)| 2 + |V’b(i’)| 2 ] is the scalar pairing function, the effective 
polarized 4-vector potential in (2+1) dimensions (so reduced to 3 components) is 
A^(r) = (U/q) [|?Mr)| 2 - | V +( r )| 2 - n/U, |^(r)| 2 - |VtB(r)| 2 , 0], and q is an effective 
charge. As before 1 2 is the two-dimensional unit matrix. Equations (107) comprise 


the relativistic Bogoliubov-de Gennes equations also known as the Dirac-Bogoliubov-de 
Gennes equations gZl 08]. In the special case of a uniform condensate which solves 


the nonlinear Dirac equation we have |Vki| 2 = \^b \ 2 = h/U and Eq. (107) yields the 
eigenvalues 


Eh — A 1 


yj ( hcik ) 2 + ( m s c 2 ) 2 ± ( m s c 2 + p) 


+ 4/i 2 


(108) 


where the magnitude of the quasiparticle momentum fc = |k| labels the eigenstates. The 
signs outside of the radical relate to pseudospin valley states and those inside the radical 
to the particle-hole Narnbu states. The spectrum Eq. (108) is plotted in Fig. [9} 



Figure 9. (color online) Relativistic BCS spectrum. Four branches corresponding to 


the sign combinations in Eq. (108): ++, -|—,-, —K (top to bottom). The vertical 


scale is in units of the chemical potential and the horizontal scale is in units of the 
reciprocal of the condensate healing length. We have indicated the superconducting 
gap 2A p located at k gap —(iJ./tici) [l + (2 m s cf /g) 2 ] 1 ^ 2 . 


In our BCS analogy the electron-positron spectral gap is 2 m s c 2 and the 
superconducting gap is 2A p = Ap located at k 2 ap =kp[l + 4(q/t+) 2 ]. Here the analogs of 
the Fermi wavenumber and velocity are = p/hci and Vf = hkF/m s = p/m s Ci. In our 
analogy we see that the Fermi momentum is inversely related to the relativistic healing 
length £ = hci/p (defined in terms of chemical potential p) such that pp = h/£. Consider 
the non-superconducting limit of the spectrum in which A p —> 0 =$> p/U <g. 1, keeping 


in mind that in the analog BCS system the Cooper pair mass is A p /c 2 . Equation (108) 
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Ek 


± 


\J ( hcik ) 2 + (m s c 2 ) 2 ± 


m s c l 


(109) 


Furthermore, for excitations much smaller than the mass gap, hcik -C m s cf, the four 
branches of the spectrum reduce to two free positive and negative energy Schrodinger- 
like excitations 


Ek 


2m s ci J 


( 110 ) 


and two similar excitations but shifted by constant potentials ±2m s c 2 

2 


Eu = ± 


he, 


2 m s c, 


k 2 + 2 m s c l 


( 111 ) 


Conversely, when hc,k m s cj linear propagation dominates the spectrum in which case 
Eq. (109) reduces to 

E k = ± ( hc,k ± m s c 2 ) , (112) 


corresponding to two copies of the Dirac spectrum shifted up or down by ±m s c 2 . 
Conversely, in the strongly superconducting regime the pairing function, and hence 
Cooper pair mass, is large compared to the positron-electron mass. This condition 
reads A p m s c 2 or 2p m s cf, in terms of the chemical potential of the condensate. 
In the limit where the kinetic energy is large as well, i.e., hcik ^ m s c 2 , the four branches 
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Figure 10. (color online) Schematic of regimes for the augmented RLSE. (a) Limiting 
theories for relativistic and non-relativistic regimes (vertical axis) versus the relative 
strengths of on-site and nearest-neighbor interactions (horizontal axis). The size of 
the superconducting gap is the same as that of the electron-positron spectral gap. 
(b) Limiting theories for relativistic and non-relativistic regimes (vertical direction) 
versus the relative strengths of superconducting and electron-positron gaps (horizontal 
direction). Here the strengths of the on-site and nearest-neighbor interactions are 
equal. In both (a) and (b) the dashed lines indicate separation between the different 
regimes and correspond to zero-temperature crossovers in the underlying lattice theory. 
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(113) 


(114) 


In the non-relativistic limit d/q = hk/m s Ci <C 1, Eq. (108) reduces to two non¬ 
propagating modes 


E k = ±yj (2 m s c\ + p) 2 + 4/i 2 , (115) 

and two modes that correspond to the particle and holes states of standard BCS theory 


E k — E i 


h 2 k 2 

2m, 


- h 


A 2 

4-1 p ■ 


(116) 


where we have reinserted the superconducting gap notation A p = 2 fi. The various limits 


and regimes are displayed in Fig. 10 


6.3. Realization in spin-orbit coupled Bose-Einstein condensates 

Our results so far can be implemented for a spin-orbit coupled BEC by considering a 2D 
pseudospin-1/2 Rashba system with variable pseudospin interactions [SB]- We note that 
the spin-orbit coupling of the Rashba form has yet to be realized but many proposed 
methods exist (see for example Refs. m Ea m eqi [9T] ). The defining Hamiltonian 
H = Hn + reads 


H 0 = [d 2 r (p 2 + 2k, p • a + hSa z ) T , 

J 2m 

Hint = (si n\ + 92 n% + 2g u ftt n 2 ) , 


(117) 

(118) 


where T = (Ti, T 2 ) T , a = {a x , a y }, h 1 ( 2 ) = |'f'i( 2 )| 2 ! d denotes the laser detuning from 
Raman resonance, and g \, g 2 , r/ 12 , are the couplings between pseudospin components. 
Note that the strength of spin-orbit coupling k depends on the relative incident angle of 
the Raman beams. In the present context we use the standard pseudospin notation 
which maps to the honeycomb lattice notation by 4 >i( 2 ) —> Ta(s). In practice, 
there are several ways to eliminate the quadratic dispersion in Eq. (117). The most 


straightforward approach would be to consider symmetric wavepackets with (p) = 0 
and momentum width A p -C 2 hn, /Mm, in which case one may safely neglect the p 2 
term in H 0 [9(2]. A second approach is to implement a setup similar to that in Ref. [ 93] . 
In this method, atoms are pumped from the two ground states via a complex external 
potential into the p = 0 state. The key point here is to maintain a stable population 
inversion given that the actual ground state is centered on a finite value of p. Both 
approaches effectively convert H to a CVE -symmetric Hamiltonian. 
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The underlying map that connects the spin-orbit coupled Hamiltonian to the 
linearized gapless RLSE requires setting the detuning to 5 = 0 in Eqs. (llT)-( 118), 
the ratio of spin-orbit coupling strength to atomic mass equal to the effective speed 
of light in the lattice k/tti = c\, and the couplings g\ = g 2 = U > 0, g \2 = 0. The 
mathematical steps of Sec. [4] may then be implemented in the same way to arrive 
at Eq. (77). In contrast, the reductions in Sec. 6.1 then require a finite detuning 
set to 5 = t 0 . Similarly, the Dirac-Bogoliubov-de-Gennes equations, Eq. (107), are 
derived from the spin-orbit coupled Hamiltonian by retaining a finite detuning and 
setting the values for the couplings g \, g 2 , and g± 2 that correctly reproduce the different 
Lorentz invariant interactions. We illustrate this for the scalar-scalar and vector-vector 
interactions, i.e., Gross-Neveu [93] and Thirring [95] models, respectively. The scalar- 


scalar interaction reads H mt = /d 2 r D(TT ) 2 = f d 2 r U (|Ti | 4 + |'I f 2 | 4 — 2|Ti|-|T 2 | 2 ), 
which is obtained in Eq. (118) for g\ — g 2 = U , g\ 2 = —2 U, where U > 0. The vector- 
vector case reads H int = Jd 2 r ^/(Ty^'h ) 2 = fd 2 rU (IT!! 4 -!- |H7 2 | 4 + 6 |\I / 1 |“| X I / 2 | 2 ), which 
requires purely repulsive interactions g\ — g 2 = U and g\ 2 = 6U. Once the particular 
form of the nonlinearity is constructed, either by adjusting on-site and nearest neighbor 
interactions for the lattice or by tuning the pseudospin interactions for the case of spin- 
orbit coupling, one may invoke the full results of Sec. [| through Secs. |6.1| and | 6.2 


In summary, our analysis in this paper describes in detail the structure of low-energy 
fluctuations near a Dirac point of a honeycomb lattice or the zero-momentum point of 
a spin-orbit coupled BEC. 


7. Conclusion 

In this article we have delineated the various constraints required for stabilizing a 
BEC at Dirac points of a honeycomb optical lattice. Energetically, we find that 
the Bose gas must be weakly interacting with excitations in the transverse direction 
suppressed relative to longitudinal ones. The latter condition can be implemented 
by using a relatively small vertical trap size. Additionally, Bloch states for the Bose 
gas must remain near enough to the Dirac point crystal momentum so that second- 
order band distortions are negligible. This condition is equivalent to the requirement 
that quasiparticle momenta remain much less than the Dirac point crystal momentum. 
Length constraints include a large quasi-2D effective healing length relative to the 
lattice spacing so that a continuum theory is physically sensible. Atomic and lattice 
parameters are related primarily by imposing the usual Landau criterion for dynamical 
stability, which relates the effective speed of light (lattice parameters) to the quasi-2D 
renormalized speed of sound (atomic parameters). 

We performed a detailed analysis of lifetimes for nonlinear Dirac vortices, 
elucidating the low-energy landscape for each solution type. Vortex lifetimes 
were computed based on dynamical instabilities induced by quantum fluctuations: 
complex eigenvalues appear in the linear spectrum for all vortex types. These 
include a complex topological vortex, topological vortices with generic winding, ring- 
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vortex, ring-vortex/soliton, vortex/soliton, Mermin-Ho, Anderson-Toulouse, and half- 
quantum vortices. The longest lived vortices are the ring-vortex, ring-vortex/soliton, 
vortex/soliton, and Anderson-Toulouse vortex with lifetimes 0.5295 s, 4.043 s, 3.841 s, 
and 4.041 s, respectively. 

A significant part of our work was devoted to the derivation and analysis of the 
relativistic linear stability equations (RLSE). We demonstrated that the RLSE reduce 
to several well known equations. The presence of a mass gap through an offset in 
the sublattice potential depths allows for an interpolation between the RLSE and 
the BdGE. By tuning the ratio of the gap to the chemical potential between small 
and large values, the governing equations for quasiparticles vary continuously between 
RLSE and Bogoliubov-de Gennes equations (BdGE) passing through a Klein-Gordon 
type structure associated with fluctuations of the nonlinear Klein-Gordon equation. In 
the particle regime where momenta are large compared to the interaction strength, 
the three types of stability equations reduce to the standard Dirac and Schrodinger 
equations with the Klein-Gordon equation interpolating between these. In the single¬ 
mode approximation, where the pseudospin valley spatial functions are equal, the RLSE 
reduce to the Andreev equations for electrons in inhomogeneous superconductors. For 
zero-energy modes residing at the core of a defect such as a vortex, the RLSE reduce to 
the Majorana equation with the Majorana mass determined by the local density of the 
condensate at the “notch” in the case of the vortex/soliton, and equal to the chemical 
potential in the general case of higher winding vortices (£ > 1). 

By including nearest-neighbor interactions and a mass gap we have shown that the 
RLSE transform to the Dirac-Bogoliubov-de-Gennes equations, which describe Cooper 
pairing of relativistic fermions. The additional Nambu space elevates the two-spinors in 
two spatial dimensions to a four component object consistent with our RLSE. The non- 
relativistic limit is defined for quasiparticle momenta much smaller than the momentum 
scale set by the mass gap, in which case we recover standard BCS theory. In the analog 
picture the BCS pairing function is mapped to the total local condensate density, that 
is, the sum of squared moduli of the sublattice amplitudes. Superconductivity is strong 
or weak depending on the magnitude of the pairing function relative to the mass gap 
energy. We have shown that when the pairing function transforms as a scalar under 
the Lorentz group the absence of internal structure for the scalar term leaves an extra 
degree of freedom in the form of a vector potential. The difference in sublattice densities 
acts as an additional polarized vector potential acting on the pseudospin-Nambu spinor. 

Interesting research directions that extend the work presented in this article could 
include elevating the boson-honeycomb lattice problem to a relativistic field theory. The 
lowest-band approximation would still be viable provided the theory is regularized by 
imposing an upper momentum cutoff at the lattice scale. The various classes of Lorentz 
quartic interactions may be constructed by including nearest-neighbor interactions in 
the lattice, as we have outlined in Sec. [6] of this article. It has been demonstrated 
that quartic interactions are fundamentally constrained by the conformal structure of 
all the terms of a particular relativistic Lagrangian [96]. Thus, by tuning the sign and 
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strength of nearest-neighbor interactions it may be possible to observe quantum phase 
transitions in the superfluid phase between different conformal theories associated with 
various relativistic field theories. 
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